Author

Diego Pinto 

── Attaching core tidyverse packages ──────────────────────── tidyverse 2.0.0 ──
✔ dplyr     1.1.4     ✔ readr     2.1.5
✔ forcats   1.0.0     ✔ stringr   1.5.1
✔ lubridate 1.9.3     ✔ tibble    3.2.1
✔ purrr     1.0.2     ✔ tidyr     1.3.1
── Conflicts ────────────────────────────────────────── tidyverse_conflicts() ──
✖ dplyr::filter() masks stats::filter()
✖ dplyr::lag()    masks stats::lag()
ℹ Use the conflicted package (<http://conflicted.r-lib.org/>) to force all conflicts to become errors
Loading required package: viridisLite


Attaching package: 'cowplot'


The following object is masked from 'package:lubridate':

    stamp



Attaching package: 'gridExtra'


The following object is masked from 'package:dplyr':

    combine



Attaching package: 'psych'


The following objects are masked from 'package:ggplot2':

    %+%, alpha
Rows: 503 Columns: 140
── Column specification ────────────────────────────────────────────────────────
Delimiter: ","
chr   (8): school, race, hispanic_latino, gender, currentSchoolMusic, school...
dbl (132): survey_id, grade, s1_1, s1_2, s1_3, s1_4, s1_5r, s1_6, s1_7, s1_8...

ℹ Use `spec()` to retrieve the full column specification for this data.
ℹ Specify the column types or set `show_col_types = FALSE` to quiet this message.

Density Plots

Code 
ggplot(data=data, aes(x=msc_full, 
                      group=currentSchoolMusic, 
                      fill=currentSchoolMusic)) +
    geom_density(adjust=1.5, alpha=.4) +
    theme_ipsum()  +
    ggtitle("Enrollment by Elective Class")

Code 
ggplot(data=data, aes(x=msc_full, 
                      group=school_music_elective, 
                      fill=school_music_elective)) +
    geom_density(adjust=1.5, alpha=.4) +
    theme_ipsum()  +
    ggtitle("Global Music Self-Concept by Elective Category")

Code 
ggplot(data=data, aes(x=musical_ability_msc, 
                      group=school_music_elective, 
                      fill=school_music_elective)) +
    geom_density(adjust=1.5, alpha=.4) +
    theme_ipsum()  +
    ggtitle("Academic Music Self-Concept by Elective Category")

Code 
ggplot(data=data, aes(x=msc_no_ability, 
                      group=school_music_elective, 
                      fill=school_music_elective)) +
    geom_density(adjust=1.5, alpha=.4) +
    theme_ipsum()  +
    ggtitle("Non-Academic Music Self-Concept by Elective Category")

Code 
ggplot(data=data, aes(x=msc_full, 
                      group=currentSchoolMusic, 
                      fill=currentSchoolMusic)) +
    geom_density(adjust=1.5) +
    theme_ipsum() +
    facet_wrap(~currentSchoolMusic) +
    theme(legend.position="none",
          panel.spacing = unit(0.1, "lines"),
          axis.ticks.x=element_blank()
    ) +
  ggtitle ("Enrollment by Elective Class")

Code 
mindset_plot_with_legend <- ggplot(data = data, aes(x = growth_mindset, 
                        group = school_music_elective, 
                        fill = school_music_elective)) +
  geom_density(adjust = 1.5, alpha = .8) +
  scale_fill_manual(values = c("No" = "#4E2A84", "Yes" = "#E4E0EE")) +
  theme_ipsum() +
  labs(x = "Mindset of Music Ability", fill = "Enrolled in School Music?") +
  theme(legend.position = "bottom", legend.margin = margin(l = -10, unit = "pt"))
Code 
mindset_plot <- ggplot(data = data, aes(x = growth_mindset, 
                        group = school_music_elective, 
                        fill = school_music_elective)) +
  geom_density(alpha = .8) + #removed "adjust = 1.5" to show more nuance on the density line
  scale_fill_manual(values = c("No" = "#4E2A84", "Yes" = "#E4E0EE")) +
  theme_ipsum() +
  labs(x = "Mindset of Music Ability", fill = "Enrolled in School Music?")  +
  theme(legend.position = "none", plot.title = element_text(size = 10))
Code 
active_engagement_plot <- ggplot(data = data, aes(x = active_engagement, 
                        group = school_music_elective, 
                        fill = school_music_elective)) +
  geom_density(alpha = .8) +  #removed "adjust = 1.5" to show more nuance on the density line
  scale_fill_manual(values = c("No" = "#4E2A84", "Yes" = "#E4E0EE")) +
  theme_ipsum() +
  labs(x = "Active Engagement", fill = "Enrolled in School Music?") +
  theme(legend.position = "none", plot.title = element_text(size = 10))
Code 
#Combine plots
combined_plot <- grid.arrange(mindset_plot, active_engagement_plot, ncol = 2)

Code 
# plot1 with modified legend
plot1_legend <- ggplot(data, 
                       aes(x = growth_mindset, y = school_music_elective, col = school_music_elective)) + 
  geom_point(size = 4) + 
  scale_color_manual(values = c("No" = "#4E2A84", "Yes" = "#E4E0EE"), name = "Enrolled in School Music?") +
  theme(legend.position = "bottom")

# function to extract legend from plot
get_only_legend <- function(plot) { 
  plot_table <- ggplot_gtable(ggplot_build(plot))
  legend_plot <- which(sapply(plot_table$grobs, function(x) x$name) == "guide-box")
  legend <- plot_table$grobs[[legend_plot]]
  return(legend)
}

# extracting modified legend from plot1 using the above function
legend <- get_only_legend(plot1_legend)

# final combined plot with shared modified legend
grid.arrange(combined_plot, legend, nrow = 2, heights = c(10, 1))

Histogram

Code 
data %>% 
  ggplot( aes(x=msc_full)) +
    geom_histogram(binwidth = 4, fill = "#4E2A84", color = "#B6ACD1", alpha= 0.9) +
  ggtitle("Global Music Self-Concept") +
    theme_ipsum()  +
    theme(plot.title = element_text(size=15)
    )

Code 
#Distribution of Music Self-Concept Subscales and Mindset
data %>% 
  select(active_engagement_z:msc_full_z) %>% 
  multi.hist() +
    theme_ipsum()

NULL
Code 
#Distribution of Music Participation Variables
data %>% 
  select(schoolChoir:selfTaught) %>% 
  multi.hist() +
    theme_ipsum()

NULL
Code 
filtered_data <- subset(data, school_music_elective == "Yes")

ggplot(filtered_data, aes(x = growth_mindset)) +
  geom_histogram(binwidth = 1, fill = "#4E2A84", color = "white", alpha = 0.7) +
  theme_minimal() +
  labs(title = "Mindset of Students Enrolled in Elective Music Classes",
       x = "Mindset Score",
       y = "Frequency")

Box Plot

Code 
data %>%
  ggplot( aes(x=currentSchoolMusic, y=msc_full, fill=currentSchoolMusic)) +
    geom_boxplot() +
    scale_fill_viridis(discrete = TRUE, alpha=0.6) +
    geom_jitter(color="black", size=0.4, alpha=0.9) +
    theme_ipsum() +
    theme(
      legend.position="none",
      plot.title = element_text(size=11)
    ) +
    ggtitle("Distribution by Elective Class") +
    xlab("")

Violin Plots

Code 
data %>%
  ggplot( aes(x=currentSchoolMusic, y=msc_full, fill=currentSchoolMusic)) +
    geom_violin() +
    scale_fill_viridis(discrete = TRUE, alpha=0.6) +
    geom_jitter(color="black", size=0.4, alpha=0.9) +
    theme_ipsum() +
    theme(
      legend.position="none",
      plot.title = element_text(size=11)
    ) +
    ggtitle("Global Music Self-Concept per Elective Class") +
    xlab("")

Code 
data %>%
  ggplot( aes(x=currentSchoolMusic, y=musical_ability_msc, fill=currentSchoolMusic)) +
    geom_violin() +
    scale_fill_viridis(discrete = TRUE, alpha=0.6) +
    geom_jitter(color="black", size=0.4, alpha=0.9) +
    theme_ipsum() +
    theme(
      legend.position="none",
      plot.title = element_text(size=11)
    ) +
    ggtitle("Academic Music Self-Concept per Elective Class") +
    xlab("")

Code 
data %>%
  ggplot( aes(x=currentSchoolMusic, y=msc_no_ability, fill=currentSchoolMusic)) +
    geom_violin() +
    scale_fill_viridis(discrete = TRUE, alpha=0.6) +
    geom_jitter(color="black", size=0.4, alpha=0.9) +
    theme_ipsum() +
    theme(
      legend.position="none",
      plot.title = element_text(size=11)
    ) +
    ggtitle("Non-Academic Music Self-Concept per Elective Class") +
    xlab("")

Radar Charts

Full Sample

Code 
data_subset <- data %>% 
  select(growth_mindset_z:movement_dance_msc_z)

average_values <- colMeans(data_subset)


max_values <- apply(data_subset, 2, max)
min_values <- apply(data_subset, 2, min)

data_for_radar <- as.data.frame(rbind(max_values, min_values, average_values))


radarchart(data_for_radar)

Code 
radarchart(data_for_radar, axistype = 1, 
           pcol = "#4E2A84",        # Color of the average values
           plwd = 2,              # Line width
           plty = 1,              # Line type
           cglcol = "grey",       # Color of the grid lines
           cglty = 1,             # Type of the grid lines
           cglwd = 0.8,           # Width of the grid lines
           axislabcol = "grey",   # Color of the axis labels
           vlcex = 0.8)           # Control the size of variable labels

Gender

Code 
# Set graphic colors
library(RColorBrewer)
coul <- brewer.pal(3, "BuPu")
colors_border <- coul
library(scales)

Attaching package: 'scales'
The following objects are masked from 'package:psych':

    alpha, rescale
The following object is masked from 'package:viridis':

    viridis_pal
The following object is masked from 'package:purrr':

    discard
The following object is masked from 'package:readr':

    col_factor
Code 
colors_in <- alpha(coul,0.3)

data_subset_gender <- data %>% 
  select(gender, growth_mindset_z:movement_dance_msc_z) %>%
  filter(gender %in% c("Male", "Female")) %>% 
  rename("Growth Mindset" = growth_mindset_z,
         "Movement & Dance" = movement_dance_msc_z,
         "Ideal Musical Self" = ideal_musical_self_msc_z,
         "Adaptive Musical Self" = adaptive_musical_self_msc_z,
         "Music Ability" = musical_ability_msc_z,
         "Community" = community_msc_z,
         "Mood Management" = mood_management_msc_z)

# Group the data by 'gender' and calculate the mean for each variable
average_by_gender <- data_subset_gender %>%
  group_by(gender) %>%
  summarise(across(everything(), mean, na.rm = TRUE))
Warning: There was 1 warning in `summarise()`.
ℹ In argument: `across(everything(), mean, na.rm = TRUE)`.
ℹ In group 1: `gender = "Female"`.
Caused by warning:
! The `...` argument of `across()` is deprecated as of dplyr 1.1.0.
Supply arguments directly to `.fns` through an anonymous function instead.

  # Previously
  across(a:b, mean, na.rm = TRUE)

  # Now
  across(a:b, \(x) mean(x, na.rm = TRUE))
Code 
# Check the average data by gender
print(average_by_gender)
# A tibble: 2 × 8
  gender `Growth Mindset` `Mood Management` Community `Music Ability`
  <chr>             <dbl>             <dbl>     <dbl>           <dbl>
1 Female           -0.379              1.37     0.519          -0.244
2 Male              0.390             -1.49    -0.507           0.316
# ℹ 3 more variables: `Adaptive Musical Self` <dbl>,
#   `Ideal Musical Self` <dbl>, `Movement & Dance` <dbl>
Code 
# Create the max and min rows based on the entire dataset
max_values <- apply(data_subset_gender[, -which(names(data_subset_gender) == "gender")], 2, max)

min_values <- apply(data_subset_gender[, -which(names(data_subset_gender) == "gender")], 2, min)

# Combine the max, min, and average values for the two gender groups into a dataframe
data_for_radar <- rbind(max_values, min_values, as.data.frame(average_by_gender[,-1]))

# Ensure the data is a dataframe and proper column names are retained
data_for_radar <- as.data.frame(data_for_radar)

# Create a radar chart with multiple groups (one for each gender)
radarchart(data_for_radar, axistype = 1, 
           pcol = c("#4E2A84", "#B6ACD1"),# Different colors for each gender
           plwd = 2,                      # Line width
           plty = 1,                      # Line type
           cglcol = "grey",               # Color of the grid lines
           cglty = 1,                     # Type of the grid lines
           cglwd = 0.8,                   # Width of the grid lines
           axislabcol = "grey",           # Color of the axis labels
           vlcex = 0.8)                   # Control the size of variable labels

# Add a legend to distinguish between genders
legend("topright", legend = unique(data_subset_gender$gender), col = c("#4E2A84", "#B6ACD1"), lty = 1, lwd = 2)

Electives

Code 
data_subset_elective <- data %>% 
  select(currentSchoolMusic, growth_mindset_z:movement_dance_msc_z) %>%
  rename("Growth Mindset" = growth_mindset_z,
         "Movement & Dance" = movement_dance_msc_z,
         "Ideal Musical Self" = ideal_musical_self_msc_z,
         "Adaptive Musical Self" = adaptive_musical_self_msc_z,
         "Music Ability" = musical_ability_msc_z,
         "Community" = community_msc_z,
         "Mood Management" = mood_management_msc_z)

# Group the data by 'currentSchoolMusic' and calculate the mean for each variable
average_by_elective <- data_subset_elective %>%
  group_by(currentSchoolMusic) %>%
  summarise(across(everything(), mean, na.rm = TRUE))

# Check the average data by elective
print(average_by_elective)
# A tibble: 4 × 8
  currentSchoolMusic `Growth Mindset` `Mood Management` Community
  <chr>                         <dbl>             <dbl>     <dbl>
1 Band                         1.64              -0.796    -0.633
2 Choir                        0.861              0.662     0.765
3 No music                    -0.912             -0.192    -0.188
4 Other                       -0.0996             1.17      0.519
# ℹ 4 more variables: `Music Ability` <dbl>, `Adaptive Musical Self` <dbl>,
#   `Ideal Musical Self` <dbl>, `Movement & Dance` <dbl>
Code 
# Create the max and min rows based on the entire dataset
max_values <- apply(data_subset_elective[, -which(names(data_subset_elective) == "currentSchoolMusic")], 2, max)

min_values <- apply(data_subset_elective[, -which(names(data_subset_elective) == "currentSchoolMusic")], 2, min)

# Combine the max, min, and average values for the gender groups into a dataframe
data_for_radar <- rbind(max_values, min_values, as.data.frame(average_by_elective[,-1]))

# Ensure the data is a dataframe and proper column names are retained
data_for_radar <- as.data.frame(data_for_radar)

# Create a radar chart with multiple groups (one for each elective)
radarchart(data_for_radar, axistype = 1, 
           pcol = c("red", "blue", "green", "purple"),# Different colors for each elective
           plwd = 2,                      # Line width
           plty = 1,                      # Line type
           cglcol = "grey",               # Color of the grid lines
           cglty = 1,                     # Type of the grid lines
           cglwd = 0.8,                   # Width of the grid lines
           axislabcol = "grey",           # Color of the axis labels
           vlcex = 0.8)                   # Control the size of variable labels

# Add a legend to distinguish between electives
legend("topright", legend = unique(data_subset_elective$currentSchoolMusic), col = c("red", "blue", "green", "purple"), lty = 1, lwd = 2)

Music Elective (Yes or No) $1

Code 
data_subset_music_elective <- data %>% 
  select(school_music_elective, growth_mindset, mood_management_msc:movement_dance_msc) %>%
  rename("Growth Mindset" = growth_mindset,
         "Movement & Dance" = movement_dance_msc,
         "Ideal Musical Self" = ideal_musical_self_msc,
         "Adaptive Musical Self" = adaptive_musical_self_msc,
         "Music Ability" = musical_ability_msc,
         "Community" = community_msc,
         "Mood Management" = mood_management_msc)

# Group the data by 'school_music_elective' and calculate the mean for each variable
average_by_music_elective <- data_subset_music_elective %>%
  group_by(school_music_elective) %>%
  summarise(across(everything(), mean, na.rm = TRUE))

# Check the average data by music_elective
print(average_by_music_elective)
# A tibble: 2 × 8
  school_music_elective `Growth Mindset` `Mood Management` Community
  <chr>                            <dbl>             <dbl>     <dbl>
1 No                                11.8              19.3      9.42
2 Yes                               13.3              19.6      9.76
# ℹ 4 more variables: `Music Ability` <dbl>, `Adaptive Musical Self` <dbl>,
#   `Ideal Musical Self` <dbl>, `Movement & Dance` <dbl>
Code 
# Create the max and min rows based on the entire dataset
max_values <- apply(data_subset_music_elective[, -which(names(data_subset_music_elective) == "school_music_elective")], 2, max)

min_values <- apply(data_subset_music_elective[, -which(names(data_subset_music_elective) == "school_music_elective")], 2, min)


data_for_radar_music_elective <- as.data.frame(bind_rows(max_values, min_values, average_by_music_elective [,-1]))


# Combine the max, min, and average values for the two music_elective groups into a dataframe
#data_for_radar_music_elective <- rbind(max_values, min_values, as.data.frame(average_by_music_elective[,-1]))

# Ensure the data is a dataframe and proper column names are retained
#data_for_radar_music_elective <- as.data.frame(data_for_radar_music_elective)

row.names(data_for_radar_music_elective) <- c("Max", "Min", unique(data_subset_music_elective$school_music_elective))


# Create a radar chart with multiple groups (one for each music_elective)
radarchart(data_for_radar_music_elective, axistype = 1, 
           pcol = c("#4E2A84", "#B6ACD1"),# Different colors for each elective
           plwd = 2,                      # Line width
           plty = 1,                      # Line type
           cglcol = "grey",               # Color of the grid lines
           cglty = 1,                     # Type of the grid lines
           cglwd = 0.8,                   # Width of the grid lines
           axislabcol = "grey",           # Color of the axis labels
           vlcex = 0.8)                   # Control the size of variable labels

# Add a legend to distinguish between genders
legend("bottomright", legend = rownames(data_for_radar_music_elective[-c(1,2),]), col = c("#4E2A84", "#B6ACD1"), lty = 1, lwd = 2, title = "School Music?")

  # Add title
title(main = "Average by Enrollment in School Music")

Music Elective (Yes or No) #2

Code 
data_subset_music_elective <- data %>% 
  select(school_music_elective, growth_mindset, mood_management_msc:movement_dance_msc) %>%
  rename("Growth Mindset" = growth_mindset,
         "Movement & Dance" = movement_dance_msc,
         "Ideal Musical Self" = ideal_musical_self_msc,
         "Adaptive Musical Self" = adaptive_musical_self_msc,
         "Music Ability" = musical_ability_msc,
         "Community" = community_msc,
         "Mood Management" = mood_management_msc)

# Group the data by 'school_music_elective' and calculate the mean for each variable
average_by_music_elective <- data_subset_music_elective %>%
  group_by(school_music_elective) %>%
  summarise(across(everything(), mean, na.rm = TRUE))

# Check the average data by music_elective
print(average_by_music_elective)
# A tibble: 2 × 8
  school_music_elective `Growth Mindset` `Mood Management` Community
  <chr>                            <dbl>             <dbl>     <dbl>
1 No                                11.8              19.3      9.42
2 Yes                               13.3              19.6      9.76
# ℹ 4 more variables: `Music Ability` <dbl>, `Adaptive Musical Self` <dbl>,
#   `Ideal Musical Self` <dbl>, `Movement & Dance` <dbl>
Code 
# Calculate max and min values
max_values <- apply(data_subset_music_elective[, -which(names(data_subset_music_elective) == "school_music_elective")], 2, max)
min_values <- apply(data_subset_music_elective[, -which(names(data_subset_music_elective) == "school_music_elective")], 2, min)

# Convert max and min values to data frames
max_values_df <- as.data.frame(t(max_values), stringsAsFactors = FALSE)
min_values_df <- as.data.frame(t(min_values), stringsAsFactors = FALSE)

# Set explicit row names for max and min
row.names(max_values_df) <- "Max"
row.names(min_values_df) <- "Min"

# Prepare average values without the first column (school_music_elective)
average_values_df <- as.data.frame(average_by_music_elective[,-1], stringsAsFactors = FALSE)

# Store the row names for the factors (Yes, No)
factor_levels <- unique(data_subset_music_elective$school_music_elective)

# Set row names for the average values
row.names(average_values_df) <- factor_levels


# Combine max, min, and average values into one data frame
data_for_radar_music_elective <- rbind(max_values_df, min_values_df, average_values_df)

# Now reassign the correct row names explicitly to avoid any swap
row.names(data_for_radar_music_elective) <- c("Max", "Min", factor_levels)

# Check the resulting data frame
print(data_for_radar_music_elective)
    Growth Mindset Mood Management Community Music Ability
Max       16.00000        24.00000 16.000000      20.00000
Min        4.00000         6.00000  4.000000       5.00000
Yes       11.81992        19.27969  9.417625      10.93870
No        13.30165        19.57025  9.760331      13.50413
    Adaptive Musical Self Ideal Musical Self Movement & Dance
Max              16.00000           20.00000        16.000000
Min               4.00000            5.00000         4.000000
Yes              11.70881           12.14559         9.808429
No               12.17355           13.14463        10.235537
Code 
# Create a radar chart with multiple groups (one for each music_elective)
radarchart(data_for_radar_music_elective, axistype = 1, 
           pcol = c("#4E2A84", "#B6ACD1"),# Different colors for each elective
           plwd = 2,                      # Line width
           plty = 1,                      # Line type
           cglcol = "grey",               # Color of the grid lines
           cglty = 1,                     # Type of the grid lines
           cglwd = 0.8,                   # Width of the grid lines
           axislabcol = "grey",           # Color of the axis labels
           vlcex = 0.8)                   # Control the size of variable labels

# Add a legend to distinguish between genders
legend("bottomright", legend = rownames(data_for_radar_music_elective[-c(1,2),]), col = c("#4E2A84", "#B6ACD1"), lty = 1, lwd = 2, title = "School Music?")

  # Add title
title(main = "Average by Enrollment in School Music")

Bar Plots

Music Class Enrollment (total and by school)

Code 
#Count and Proportion per class
music_elective_table <- data %>% 
 count(currentSchoolMusic)%>% 
  mutate("%" = round(n/sum(n)*100,2))

#Proportion per school per class
 round(prop.table(table(data$school,data$currentSchoolMusic),1)*100, 2)
   
     Band Choir No music Other
  A 26.40 31.20    18.40 24.00
  B 14.81 22.22    62.14  0.82
  C 15.56 12.59    64.44  7.41
Code 
 #Percentages
percentages <- data %>%
  group_by(school, currentSchoolMusic) %>%
  summarise(count = n()) %>%
  group_by(school) %>%
  mutate(percentage = count/sum(count)*100)
`summarise()` has grouped output by 'school'. You can override using the
`.groups` argument.
Code 
 #Frequency
 data %>% 
  group_by(school, currentSchoolMusic) %>% 
  summarise(n = n()) %>% 
  pivot_wider(names_from = currentSchoolMusic,
              values_from = n)
`summarise()` has grouped output by 'school'. You can override using the
`.groups` argument.
# A tibble: 3 × 5
# Groups:   school [3]
  school  Band Choir `No music` Other
  <chr>  <int> <int>      <int> <int>
1 A         33    39         23    30
2 B         36    54        151     2
3 C         21    17         87    10
Code 
ggplot(percentages, aes(fill=currentSchoolMusic, y=percentage, x=school)) + 
    geom_bar(position="stack", stat="identity") + 
  geom_text(aes(label = paste0(round(percentage), "%")), 
            position = position_stack(vjust = 0.5), 
            size = 3, 
            color = "grey", 
            fontface = "bold") +
    scale_fill_viridis(discrete = T) +
    theme_ipsum() +
    labs(x = "School", y = "%", fill = "School Music Electives")

Heat Map

Are the observed frequencies in elective participation by school significantly different from expected frequencies in the distribution?

Code 
# Perform Chi-square test
chisq_result <- chisq.test(table(data$school, data$currentSchoolMusic))

chisq_result

    Pearson's Chi-squared test

data:  table(data$school, data$currentSchoolMusic)
X-squared = 106.39, df = 6, p-value < 2.2e-16
Code 
# extracting adjusted residuals
adjusted_residuals <- residuals(chisq_result, type = "pearson")

# Converting adjusted residuals to a matrix format
adjusted_residuals_matrix <- matrix(adjusted_residuals, nrow = nrow(chisq_result$observed))

# Reshaping data for ggplot
library(reshape2)

Attaching package: 'reshape2'
The following object is masked from 'package:tidyr':

    smiths
Code 
adjusted_residuals_df <- melt(adjusted_residuals_matrix)

# Renaming levels in Var1 and Var2
adjusted_residuals_df$Var1 <- factor(adjusted_residuals_df$Var1, levels = c("1", "2", "3"), labels = c("A", "B", "C"))

adjusted_residuals_df$Var2 <- factor(adjusted_residuals_df$Var2, levels = c("1", "2", "3", "4"), labels = c("Band", "Choir", "No music", "Other music"))


# Heatmap of adjusted residuals
ggplot(data = adjusted_residuals_df, aes(x = Var2, y = Var1, fill = value)) +
  geom_tile() +
  scale_fill_gradient2(low = "blue", mid = "white", high = "red", midpoint = 0) +
  theme_minimal() +
  labs(x = "School Music Electives", y = "School", fill = "Adjusted Residuals")

Code 
ggplot(data = adjusted_residuals_df, aes(x = Var2, y = Var1, fill = value)) +
    geom_tile() +
    scale_fill_viridis() +
    theme_minimal() +
    labs(x = "School Music Ellectives", y = "School", fill = "Adjusted Residuals")